##############################################
## Replication R file
## “What Determines the Specificity of Subsidies?”
## International Studies Quarterly (2012)
##
## by Jong Hee Park
## University of Chicago
## jhp@uchicago.edu
##
## "Wed Oct 20 14:22:20 2010"
##############################################
require(foreign)
require(sandwich)
require(xtable)
master <- read.csv("StateAid.csv")
attach(master)

##############################################
## Two Stage Estimates
## 1. ysectorgdp
##############################################
first.lm <- lm(ysectorgdp ~ lagunemp + newimex + fdi + laggdpgr + eutrade +
               govleft + coalition + election + trend + as.factor(id),
               data=master, y=TRUE)
first0 <- lm(ysectorgdp ~ lagunemp + newimex + fdi + laggdpgr + eutrade +
               govleft + coalition + election + trend + id,
               data=master, y=TRUE)
master0 <- data.frame(model.matrix(first0))
n <- length(unique(master0$id))

betaFE <- matrix(coef(first.lm)[1:10], ,1)
u.hat <- y.mean <- res.mean <- coor.mean <- avemag.mean <- rep(NA, n)
first.res <- resid(first.lm)
for (i in 1:n){
  res.mean[i] <- mean(first.res[master0$id==i])
  y.mean[i] <- mean(first.lm$y[master0$id==i])
  x.mean <- matrix(apply(subset(master0, id==i), 2, mean)[1:10], 1, )
  u.hat[i] <- y.mean[i] - x.mean%*%betaFE - res.mean[i]
}
## Regress e with x using n = number of groups
for (i in 1:n){
  coor.mean[i] <- mean(master$stand.pca[master$id==i])
  avemag.mean[i] <- mean(log(master$avemag2[master$id==i])) 
}
coor.mean2 <- coor.mean^2
second.lm <- lm(u.hat ~ coor.mean + coor.mean2 + avemag.mean - 1)
second.res <- resid(second.lm)

adjust.sector.fe <- rep(as.list(NA), n)
for (i in 1:n){
  ni <- sum(master$id == i)
  if (sum(names(second.res)==i) == 0){
    adjust.sector.fe[[i]] <- rep(NA, ni)
  }
  else{
    adjust.sector.fe[[i]] <- rep(second.res[which(names(second.res)==i)], ni)
  }
}
master$adjust.sector.fe <- unlist(adjust.sector.fe)

## Use second.res as country dummies
third.lm1 <- lm(ysectorgdp ~  lagunemp + newimex + fdi +
                laggdpgr + eutrade + govleft + coalition + election + trend + log(avemag2) +
                stand.pca + stand.pca2 + adjust.sector.fe,
                data=master, y=TRUE)
third.lm2 <- lm(ysectorgdp ~  lagunemp + newimex + fdi +
                laggdpgr + eutrade + coalition + election + trend + log(avemag2) +
                govleft *(stand.pca + stand.pca2) + adjust.sector.fe,
                data=master, y=TRUE)
third.lm3 <- lm(ysectorgdp ~  lagunemp + newimex + fdi +
                laggdpgr + eutrade + govleft + coalition + election + trend + 
                log(avemag2)*(stand.pca + stand.pca2) + adjust.sector.fe,
                data=master, y=TRUE)
output1 <- cbind(summary(third.lm1)$coef[,1], sqrt(diag(vcovHAC(third.lm1))))
output2 <- cbind(summary(third.lm2)$coef[,1], sqrt(diag(vcovHAC(third.lm2))))
output3 <- cbind(summary(third.lm3)$coef[,1], sqrt(diag(vcovHAC(third.lm3))))
#
#################################
## ygeneralgdp
#################################
## Compute the residual for all covariates but x
## e = f(x, fixed effects)
first.lm <- lm(ygeneralgdp ~ lagunemp + newimex + fdi + laggdpgr + eutrade +
               govleft + coalition + election + trend + as.factor(id),
               data=master, y=TRUE)
first0 <- lm(ygeneralgdp ~ lagunemp + newimex + fdi + laggdpgr + eutrade +
               govleft + coalition + election + trend + id,
               data=master, y=TRUE)
master0 <- data.frame(model.matrix(first0))
n <- length(unique(master0$id))

betaFE <- matrix(coef(first.lm)[1:10], ,1)
u.hat <- y.mean <- res.mean <- coor.mean <- avemag.mean <- rep(NA, n)
first.res <- resid(first.lm)
for (i in 1:n){
  res.mean[i] <- mean(first.res[master0$id==i])
  y.mean[i] <- mean(first.lm$y[master0$id==i])
  x.mean <- matrix(apply(subset(master0, id==i), 2, mean)[1:10], 1, )
  u.hat[i] <- y.mean[i] - x.mean%*%betaFE - res.mean[i]
}
## Regress e with x using n = number of groups
for (i in 1:n){
  coor.mean[i] <- mean(master$stand.pca[master$id==i])
  avemag.mean[i] <- mean(log(master$avemag2[master$id==i])) 
}
coor.mean2 <- coor.mean^2
second.lm <- lm(u.hat ~ coor.mean + coor.mean2 + avemag.mean -1)
second.res <- resid(second.lm)

adjust.general.fe <- rep(as.list(NA), n)
for (i in 1:n){
  ni <- sum(master$id == i)
  if (sum(names(second.res)==i) == 0){
    adjust.general.fe[[i]] <- rep(NA, ni)
  }
  else{
    adjust.general.fe[[i]] <- rep(second.res[which(names(second.res)==i)], ni)
  }
}
master$adjust.general.fe <- unlist(adjust.general.fe)

## Use second.res as country dummies
third.lm4 <- lm(ygeneralgdp ~  lagunemp + newimex + fdi +
                laggdpgr + eutrade + govleft + coalition + election + trend + log(avemag2) +
                stand.pca + stand.pca2 + adjust.general.fe,
                data=master, y=TRUE) 
third.lm5 <- lm(ygeneralgdp ~  lagunemp + newimex + fdi +
                laggdpgr + eutrade + coalition + election + trend + log(avemag2) +
                govleft*(stand.pca + stand.pca2) + adjust.general.fe,
                data=master, y=TRUE)
third.lm6 <- lm(ygeneralgdp ~  lagunemp + newimex + fdi +
                  laggdpgr + eutrade + govleft + coalition + election + trend + 
                  log(avemag2)*(stand.pca + stand.pca2) + adjust.general.fe,
                  data=master, y=TRUE)
output4 <- cbind(summary(third.lm4)$coef[,1], sqrt(diag(vcovHAC(third.lm4))))
output5 <- cbind(summary(third.lm5)$coef[,1], sqrt(diag(vcovHAC(third.lm5))))
output6 <- cbind(summary(third.lm6)$coef[,1], sqrt(diag(vcovHAC(third.lm6))))

xtable(output1, digits=3)
xtable(output2, digits=3)
xtable(output3, digits=3)
xtable(output4, digits=3)
xtable(output5, digits=3)
xtable(output6, digits=3)

##########################################
## Excluding the Netherlands
##########################################
third.lm7 <- lm(ygeneralgdp ~  lagunemp + newimex + fdi +
                laggdpgr + eutrade + govleft + coalition + election + trend + log(avemag2) +
                 stand.pca + stand.pca2 + adjust.general.fe,
                data=subset(master, country!="Netherlands"))
third.lm8 <- lm(ygeneralgdp ~  lagunemp + newimex + fdi +
                laggdpgr + eutrade + coalition + election + trend + log(avemag2) +
                govleft*(stand.pca + stand.pca2) + adjust.general.fe,
                data=subset(master, country!="Netherlands"))
third.lm9 <- lm(ygeneralgdp ~  lagunemp + newimex + fdi +
                laggdpgr + eutrade + govleft + coalition + election + trend + 
                log(avemag2)*(stand.pca + stand.pca2) + adjust.general.fe,
                data=subset(master, country!="Netherlands"))
output7 <- cbind(summary(third.lm7)$coef[,1], sqrt(diag(vcovHAC(third.lm7))))
output8 <- cbind(summary(third.lm8)$coef[,1], sqrt(diag(vcovHAC(third.lm8))))
output9 <- cbind(summary(third.lm9)$coef[,1], sqrt(diag(vcovHAC(third.lm9))))

xtable(output7, digits=3)
xtable(output8, digits=3)
xtable(output9, digits=3)

##########################################
## Excluding the Portugal
##########################################
## Use second.res as country dummies
third.lm10 <- lm(ysectorgdp ~  lagunemp + newimex + fdi +
                laggdpgr + eutrade + govleft + coalition + election + trend + log(avemag2) +
                stand.pca + stand.pca2 + adjust.sector.fe,
                data=subset(master, country!="Portugal"))
third.lm11 <- lm(ysectorgdp ~  lagunemp + newimex + fdi +
                laggdpgr + eutrade + coalition + election + trend + log(avemag2) +
                govleft *(stand.pca + stand.pca2) + adjust.sector.fe,
                data=subset(master, country!="Portugal"))
third.lm12 <- lm(ygeneralgdp ~  lagunemp + newimex + fdi +
                laggdpgr + eutrade + govleft + coalition + election + trend + log(avemag2) +
                 stand.pca + stand.pca2 + adjust.general.fe,
                data=subset(master, country!="Portugal"))
third.lm13 <- lm(ygeneralgdp ~  lagunemp + newimex + fdi +
                laggdpgr + eutrade + coalition + election + trend + log(avemag2) +
                govleft*(stand.pca + stand.pca2) + adjust.general.fe,
                data=subset(master, country!="Portugal"))
output10 <- cbind(summary(third.lm10)$coef[,1], sqrt(diag(vcovHAC(third.lm10))))
output11 <- cbind(summary(third.lm11)$coef[,1], sqrt(diag(vcovHAC(third.lm11))))
output12 <- cbind(summary(third.lm12)$coef[,1], sqrt(diag(vcovHAC(third.lm12))))
output13 <- cbind(summary(third.lm13)$coef[,1], sqrt(diag(vcovHAC(third.lm13))))

xtable(output10, digits=3)
xtable(output11, digits=3)
xtable(output12, digits=3)
xtable(output13, digits=3)

